Broadcast FM Noodling

Now to look at some wide FM in the broadcast band. This should be the same as before, but we need to be more aggressive about filtering out the extra cruft above the pilot tone at 19kHz.

In [17]: 
%pylab inline

from gnuradio import blocks
from gnuradio import gr
from gnuradio import uhd
from gnuradio import audio
from gnuradio.filter import firdes

import scipy
import scipy.signal

Populating the interactive namespace from numpy and matplotlib
In [18]: 
F_station = int(104.5e6)
Fs=int(1e6)
gain=10
N = int(5e6)

Fc = F_station - 250000 # We capture to go a bit low of the center
In [19]: 
# Set up our flow graph
# Cribbing off http://gnuradio.org/notebooks/GNURadioforSimulation.html
#
tb = gr.top_block()

# Set up the USRP for input     
usrp_source = uhd.usrp_source(",", uhd.stream_args(cpu_format="fc32", channels=[0]))
usrp_source.set_samp_rate(Fs)
usrp_source.set_center_freq(Fc, 0)
usrp_source.set_gain(gain, 0)

# Let's try to flush out the first bunch of samples
skip_head = blocks.skiphead(gr.sizeof_gr_complex, 1000)

# Limit ourselves to N samples
head = blocks.head(gr.sizeof_gr_complex, N)

# And a sink to dump them into
sink = blocks.vector_sink_c()


tb.connect(usrp_source, skip_head, head, sink) # Can use the handy serial connect method here
In [20]: 
tb.run()
tb.stop()
usrp_source.stop()
x = np.array(sink.data())
psd(x, Fs=Fs); None
In [21]: 
specgram(xp[0:Fs], Fs=Fs); None
In [22]: 
def extract_channel(f_sample, f_center, f_bw, x):
    """Extract a channel of width f_bw at f_center, from f_sample
    
    Returns Fs,y, the new sample rate and an array of complex samples
    """
    
    # Cribbing from the GNURadio xlate, we use their block diagram
    
    my_state = {} # cheating backchannel for debug
    
    # Create a LPF for our target bandwidth
    n_taps = 100 # XXX Total random shot in the dark    
    r_cutoff = float(f_bw)/f_sample    
    base_lpf = scipy.signal.remez(n_taps, [0, r_cutoff, r_cutoff+(0.5-r_cutoff)/4, 0.5], [1,0])
    
    # Shift this LPF up to our center frequency
    fwT0 = 2.0*np.pi*(f_center-f_bw/2)/f_sample
    
    lpf = base_lpf * np.exp(1.0j*fwT0 * np.arange(0, n_taps))
    
    dec_rate = int(f_sample / f_bw)
    
    x_filtered = scipy.signal.lfilter(lpf, 1.0, x)
    
    dec_is = np.arange(0, len(x_filtered), dec_rate)
    y = x_filtered[dec_is]
    y *= np.exp(-1.0j*fwT0*dec_is)
    
    my_state["n_taps"] = n_taps
    my_state["r_cutoff"] = r_cutoff
    my_state["lpf"] = lpf
    my_state["base_lpf"] = base_lpf
    my_state["fwT0"] = fwT0
    my_state["dec_rate"] = dec_rate
    # my_state["dec_is"] = dec_is
    # my_state["x_filtered"] = x_filtered
    
    return f_sample/dec_rate, y, my_state

def decode_quad(x, gain):
    xp = x[1:] * np.conj(x[:-1])
    retval = gain * np.arctan2(np.imag(xp), np.real(xp))
    #lpf = scipy.signal.remez(30, [0.05, 0.25, 0.27, 0.5], [1,0])
    #retval = scipy.signal.lfilter(lpf, 1.0, retval)
    return retval
In [23]: 
Fs_new, y, _ = extract_channel(Fs, 350000, 200000, x)
specgram(y[0:Fs], Fs=Fs_new)
Fs_new
Out[23]:
200000

And now, to decode the actual audio from the band-limited FM signal. We can just use decode_quad for this.

In [24]: 
w = decode_quad(y, 1.0)
In [25]: 
plot(w[0:1000])
Out[25]:
[<matplotlib.lines.Line2D at 0x7f2a1cde3cd0>]
In [26]: 
figure(figsize=(14,6))
specgram(w[0:30000], Fs=Fs_new); None

We can see the audio down below 20kHz, with a strong pilot tone at 19kHz. We should be using that tone to drive a frequency locking loop of some sort, but that's too much work at the moment.

There's another tone at 2×19kHz (just below 40kHz), which is marking off where stereo information goes. And, if you look closely, you'll see another little row of signal at 3×19kHz (just below 60kHz), which is the RDS data signal. More on that later.

For now, let's low-pass the audio and filter the crap out of it to try and get rid of that pilot tone.

In [27]: 
lpf_audio = scipy.signal.remez(120, [0,15000, 16000, Fs_new/2], [1,0], Hz=Fs_new)

Fs_audio = 40000

w_audio = scipy.signal.lfilter(lpf_audio, 1.0, w)
w_audio = w_audio[range(0, len(w_audio), Fs_new/Fs_audio)]
specgram(w_audio, Fs=Fs_audio)
Fs_audio
Out[27]:
40000
In [28]: 
scale_factor = 5000 / np.max(np.abs(w_audio))

(w_audio*scale_factor).astype("int16").tofile("wbfm.raw")
In [29]: 
from subprocess import call
cmd_line = "aplay -r %d -f S16_LE -t raw -c 1 wbfm.raw" % (Fs_audio)
call(cmd_line.split(" "))
Out[29]:
0